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ABSTRACT 

The predicted mass function of dark matter halos is essential in connecting observed galaxy cluster counts 
and models of galaxy clustering to the properties of the primordial density field. We determine the mass 
function in the concordance ACDM cosmology, as well as its uncertainty, using sixteen 1024 3 -particle nested- 
volume dark-matter simulations, spanning a mass range of over five orders of magnitude. Using the nested 
volumes and single-halo tests, we find and correct for a systematic error in the friends-of-friends halo-finding 
algorithm. We find a fitting form and full error covariance for the mass function that successfully describes 
the simulations' mass function and is well-behaved outside the simulations' resolutions. Estimated forecasts 
of uncertainty in cosmological parameters from future cluster count surveys have negligible contribution from 
remaining statistical uncertainties in the central cosmology multiplicity function. There exists a potentially 
non-negligible cosmological dependence (non-universality) of the halo multiplicity function. 
Subject headings: cosmology: theory — galaxies: clusters: halos — galaxies: halos 



1. INTRODUCTION 

Collapsed, virialized dark matter halos arise from density 
peaks in the initially Gaussian primordial flu ctuation field 
( Press & Schechter 1974; Bardeen et alJll986l) . The abun- 
dance of the most massive of these halos is exponentially sen- 
sitive to the amplitude of the initial fluctuation field as well 
as the mean matter density, making observed counts of their 
abundance extremely sensitive to these properties of the den- 
sity field, as well as the d ark-energy dependen t growth rate 
of the density field (e.g., Haiman et alT l2001l) . Condensa- 
tion of gas and formation of stars within these h alos leads 
to the formation of galaxies llWhite & Reesl[l978l) . In ad- 
dition, nonlinear clustering halo-models of dark matter and 
galaxies require, as a basic component, the mass function (see 
ICoorav & Shethl2002l for a recent review). 

The analytic theory of virialized object formation 
through collapse of ov erdense regions as envisioned by 
iPress & SchechterU 19741) (hereafter PS) employs the fact that 
a uniform overdensity in the Universe will evolve as a sepa- 
rate, closed universe, initially expanding with the background, 
but then slowing and turning around to collapse and virialize. 
Since the abundance of overdensity peaks only depends on the 
fluctuation scale er, the abundance of halos can be expected to 
be universal in these units. Limitations of approximations in 
the PS model, e.g., sphericity of collapse and spatial overlap, 
led to a modi fication of the original f orm with parameters fit to 
simulations (Sheth & Tormen 1999, hereafter ST). Using the 
same simulations, Jenkins et al. (2001) (hereafter J01) aban- 
doned the form of the PS motivated mass function to better 
fit the simulations' mass range, but their functional form can- 
not be extrapolated beyond the range of the fit. J01 found the 
mass function in a to be approximately universal for several 
cosmologies at the level of ~15%. 

Here, we present a quantification of the dark matter halo 
mass function and its uncertainties with a suite of sixteen 
nested-volume 1024 3 particle dark matter simulations of the 
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concordance ACDM cosmology. We quantify the uncertainty 
and full covariance of the mass function parameters. Our 
halo-finding methodology is given in 5j2j mass function deter- 
mination and error analysis is presented in |3j along with im- 
plications for future cluster surveys' sensitivities; we present 
our conclusions in |4] 

2. numerical simulations and halo mass 
determination 

We calculate the mass function of dark matter halos aris- 
ing in a concordance ACDM model by performing numer- 
ical simulations of structure growth and halo formation. We 
us e the Hashed-Oct-Tree (H OT) algorithm, initially described 
in IWarren & Salmon! lll993l) and recently co mpared in detail 
with o ther well-known cosmology codes in iHeitmann et alJ 
(2004). W e use a per-interaction e rror bound based on the 
analysis in lSalmon & Warrenl ( 119941) . The fractional error per 
interaction is set to be no worse than 10~ 5 at a redshift of 25, 
increasing to 10~ 3 at redshifts of 5 and lower. The number 
of timesteps and Plummer smoothing lengths ranged from 
1480 steps and 2.1/r'kpc (physical) for the highest resolu- 
tion simulation, to 720 steps and 98/r'kpc (comoving) for 
the largest volume. Each simulation required about 2 x 10 17 
floating point operations, which can be computed in roughly 
60 hours on a 1024 processor parallel computer using HOT. 
Overall, the results presented here required over four exaflops 
(4 x 10 18 floating point operations). 

We model a universe with flat geometry and parameters 
p = (Sl M ,{l b ,n,h,a s ) = (0.3,0.04, 1,0.7,0.9). (1) 
Initial conditions are derived fro m the transfer funct ions as 
calculated by CMBFAST ( ISeliak & Zaldarriagal 1996b . 

In order to simultaneously reduce Poisson error, verify con- 
sistency, and resolve the widest mass-scale range available by 
our techniques, we employ nested volumes with independent 
realizations of the chosen cosmology. We simulated sixteen 
boxes of sizes 96, 135, 192, 272, 384, 543, 768, 1086, 1536, 
2172, 2583 and 3072 /z _1 Mpc, with three realizations of 384 
/r'Mpc size, and two each of 272 and 3072 /r'Mpc. After 
the aggressive requirement of a minimum of 400 particles per 
halo, we measure the mass function over five orders of mag- 
nitude in mass scale; see Fig. Q. The box size minimum is 
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FIG. 1 . — Shown are the central values of the binned mass functions from 
sixteen 1024 3 simulations of the ACDM universe as crosses, with simulations 
in different colors. The best-fit form for the mass function we calculate is in 
solid (red), the Jenkins fit dashed (purple), and the Sheth-Tormen fit in dot- 
dashed (dark grey). Goodness-of-fit is poorly judged on this extreme log 
scale; it is more clearly resolved in the linear residuals of Fig.l2l 

set by requiring the largest scale growth modes in the volume 
to remain linear, while the the maximum is set by the halo 
particle number requirement (i.e., very few massive halos in 
boxes of size > 3 n _1 Gpc have more than 400 particles). Our 
nested-volume approach allows greater mass and length-scale 
resolution than pure particle number a nd force resolution i n- 
crease within a single simulation (e.g., Springel et al. 2005). 

The friends-of-friends (FOF) method JFrenk et al J 1 1 9881) 
identifies a set of particles which are spatially associated and 
contained within a given isodensity surface defined by a link- 
ing parameter, b. The linking length is defined as ni; n k = 
fen" 1 / 3 , where n is the number density of particles. 

In practice, the number of particles per halo is not sufficient 
to suppress statistical noise in the density field represented 
by the finite number of particles. This leads to a significant 
systematic error in the estimated mass of a given halo, which 
depends on the number of particles which represent the halo. 
As an illustration of one aspect of this problem, consider an 
isothermal density distribution of total mass M represented 
by N particles within a sphere of radius R. Each particle has 
mass M p = M/N. Choosing an FOF link length h\{ n ^ defines 
an isodensity surface with value 



P(hlink) : 



(2) 



where a is "a constant of order 2" ( iFrenk et al.119 88). 

Using the FOF method on a large number of sample isother- 
mal halos with parameters M = R = 1 and nv m k = (N/1 .25)~'/ 3 
while varying N leads to the results that are presented in Ta- 
ble [2 It is clear from these results that there is a strong N- 
dependence in the determined mass of the same underlying 
density distribution. Another interesting result is that the mass 
values are converging to a value which is significantly differ- 
ent from that implied from Eq. with a = 2.0. For large N, 
the correct value for a is close to 3.1, which implies that the 
overdensity of halos being found with the commonly favored 
link parameter of b = 0.2 is about 280/%, rather than the of- 
ten qu oted value o f l&Op h llCole & Lace\ll994l:I.Tenkins et alJ 
l200ll fWhite2001. 2002). 

For a well-behaved halo finding method, one should be able 
to sub-sample a simulation and derive the same mass function 
in the regime where the statistics are robust. In the case of 
FOF, this means that one should be able to randomly pick one 



640000 0.3879 0.0014 

80000 0.3999 0.0037 

10000 0.4187 0.0092 

1250 0.4513 0.0244 



N MfOF a 

500 0.4704 0.0371 

200 0.4964 0.0555 

100 0.5222 0.0765 

50 0.5586 0.1018 



TABLE 1 

The mean and variance of FOF mass vs N for an ideal halo 



of every n particles in the simulation, run the FOF method 
with a link length y/n times as large, and find the same distri- 
bution of halos. 

However, due to the A^-dependent statistical effects de- 
scribed above, one obtains substantially different mass- 
functions from a sub-sampled distribution. While it is pos- 
sible to use the binomial distribution to determine an expres- 
sion for the A^-dependent corrections to the FOF method in the 
case of an ideal isothermal density distribution, halos found in 
simulations are complex objects which are not well-described 
by such a smooth density distribution. For this reason, we 
have chosen to implement a correction which is calibrated via 
the comparison of a simulation to the sub-sampled version of 
itself. Fortunately, the correction derived from such a proce- 
dure appears to not depend strongly on the parameters of the 
simulation. We find a correction of the form 



N com ced=N(l-N-°- 6 ), 



(3) 



for a halo with N particles to do a reasonable job of correcting 
the systematic error caused by the FOF method. 

With this correction and the understanding of the behav- 
ior of FOF for highly-resolved halos discussed above, for a 
link parameter b = 0.2 (used for all further analysis in this pa- 
per) we find that the mean density of halos with respect to 
the background in our largest volume simulation is 260 ± 50, 
while in the smallest volume this value is 340 ± 60. For this 
reason, further analysis based on the results described below 
should not assume that the overdensity of FOF halos is con- 
stant with respect to scale. The process of defining a halo 
remains the largest uncertainty in defining the mass func- 
tion. New algorithms for halo finding based on an overden- 
sity criterion rather than the use of FOF will be necessary for 
straightforward comparisons with observational data. 

3. THE MASS FUNCTION 

The relation between the multiplicity function /(er) and 
mass function n is 



/(*) : 



M dn 



(4) 



p dlna 1 ' 

where p is the mean matter density of the universe. The gen- 
eral shape of the mass function is well described by the PS 
and ST forms for the multiplicity function. We therefore use a 
similar form for the mass function as that of ST, with a disen- 
tanglement of the power-law small-mass regime from the ex- 
ponential cut-off, and removal of the arbitrary collapse scale 
S c , which was nevertheless modified by the ellipsoidal col- 
lapse scale in the ST model. The minimal-parameter multi- 
plicity function required to fit the mass-function shape mea- 
sured in the simulations is 



f(a)=A(a- a + b)e 



c/o 



(5) 



with parameters q = (A,a,b,c). This form is related to the 
multiplicit y function resu lting from the barrier shape ansatz 
of Sheth & Tormen (2002), but without an associated change 
of the exponential. 
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FIG. 2. — Shown are the residuals from the binned simulation data to the fit presented in this work as square data points of different colors per simulation. 
The Jenkins fit is the solid (purple) line, ST original fit the dashed (dark gray) line, the ST fit with parameters A,a,p free with dot-dashed line (red), and the ST 
fit with a,p free and amplitude A set to require all dark matter in halos as a triple-dot-dashed line (light gray). The binned mass function from the Virgo Hubble 
Volume simulation are the asterisk points with errors (pink). 



Using, Eq. Q, we calculate the extended maximum likeli- 
hood, A(q) for the Poisson data of the mass function fo und in 
the combined set of simulations ( Eidel man et all2 004). 



lnA(q): 
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-EE 
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(6) 



of N s simulations containing Nj bins, where n (/ is the number 
of halos in mass bin i of simulation j. In bins where n,j = 
the last term is zero. The predicted number of halos in bin ; 
of simulation j between masses m (/ i and m,^ is 

r >p ~ dn(q) 



dM, 



(7) 



where Vj is the volume of the simulation. The extended max- 
imum likelihood is appropriate for Poisson data sets such 
as the mass function, and allows for an estimation of the 
goodness-of-fit for a given model, with the number of degrees 
of freedom for the multinomial nature of the distribution given 
by N— k- 1, where N is the total number of bins and k is the 
number of fitted parameters. 

The mass function from all simulations is binned with min- 
imum logarithmic width of 0.05, with the last bin expanded 
to include a minimum of 400 halos. Higher mass halos have 
negligible statistical weight and are discarded. The low-mass 
end requires 400 particles per halo, and the FOF determined 
masses are corrected as described in ^2] 

Sample variance can be a consideration when attempting to 
use the most massive halo s available in a small volume (e.g., 
see iHu & Kravtsovll2003l) . The sample variance of halos of 
mass M in a cubic volume of side length L such as our simu- 
lations is b 2 (M)a 2 (L), where b(M) is the halo bias and <j(L) is 
the fluctuation amplitude in the cubic volume. For our simula- 
tion volume samples, the maximum sample variance is 0.2% 
of the Poisson error, which occurs for the most massive halo 
bin in our smallest volume simulation. Sample variance is 
therefore negligible here. 

We perform an analysis of the likelihood, Eq. (|6j, for the 
multiplicity function, Eq. (jSJi, to obtain the simulations' re- 
sults for the parameters q free. The best-fit of the multiplicity 
function is found with parameters 

A = 0.7234 ± 0.0043 (stat.) ± 0.003 (sys.), 
a= 1.625 ±0.019 (stat.) ±0.009 (sys.), 
b = 0.2538 ± 0.003 1 (stat.) ± 0.002 (sys.), 
c= 1.1982 ± 0.0055 (stat.) ± 0.002 (sys.), 
with corresponding la statistical errors, and estimated sys- 
tematic uncertainty. The statistical error covariance is 

^1.82 5.80 0.0730 2.04 ^ 
.. 36.7 4.35 9.85 
0.958 



C = 



0.854 
2.97 



10" 



(8) 



(9) 



Error analysis was performed with the MINUIT package of 
CERNLIB. The x 2 P er degree-of-freedom (DOF) for the best 
fit is x 2 /DOF = 524.5/429, indicating a larger scatter than 
expected from the Poisson statistics. 

The uncertainty and resulting scatter in the sampling of a 
given overdensity by the FOF halo finder likely contributes 
to the higher than expected x 2 /DOF. To assess this source 
of systematic uncertainty, we vary the FOF correction Eq. Q 
within ranges fit by single simulations, as well as alternate 
binnings of the mass function, which sample this distribution 
differently. The estimated systematic uncertainties are sub- 
dominant to statistical, though not negligible. 

The residuals of the binned mass functions from all sim- 
ulations to the expected number of halos in the given bin 
[Eq. Q] of mass-function best fit are shown in Fig. [2] In 
this figure, we also plot the binn ed mass function fr om the 
Virgo Hubble Volume simulation (Jenkins et al. 2001), which 
we derive from their particle data identically to our simula- 
tions (although the initial power spectrum for those data are 
not precisely the same as our simulations). 

For comparison, we also fit the ST form for the multi- 
plicity function. With amplitude free of the constraint that 
all dark matter be within halos, the ST parameters are A = 
0.3405 ± 0.0004, a = 0.7183 ±0.001 and p = 0.157 ±0.003, 
with y 2 / dof = 1987/430. From the derived x 2 / DOF and 
Fig. |2] modifications of ST form and the J01 form clearly 
provide a poor fit to the ACDM mass function. However, 
the ST paramet ers are similar to that fit to the bias from 
simulations of Seliak & Warren (2004) using the ST form 
as given in iMandelblmmet alJ J2004Tk Therefore, the peak- 
background split ansatz of the relation b etween halo biasing 
and the mass function may be consistent tMo & Whitdl 19961 
ISheth & Tormerll999tlShefh et alJ20'ol . Although more can 
be done to verify this statistical consistency, it is beyond the 
scope of this Letter. 

To test the claim of the multiplicity function's universal- 
ity in consistently describing varying cosmologies, we find 
the residuals for the mass function derived from the ACDM 
simulations' multiplicity function versus simulations done for 
varying cosmological parameters within ACDM. Four sepa- 
rate 1024 3 simulations were run with parameters as those in 
Eq. Q except, a 384 /z _1 Mpc box with <j 8 = 0.8; a 384 /r'Mpc 
box with Q m = 0.2; a 329 /r'Mpc box with h = 0.6 and 
cr% = 0.8; and a 1536 /r'Mpc box with Q m = 0.2 and trg = 1.2. 
Residuals are shown in Fig. [3] The models are consistent 
within the 5% level, except for the fl,„ = 0.2 and cjg = 1.2 
model, which shows departures at 20%. The departure may 
be more than a Poisson fluctuation, revealing a statistically 
significant dependence of the multiplicity function on cosmol- 
ogy. Quantifying this potential dependence is warranted, but 
beyond the scope of this work. 
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FIG. 3. — Residuals are shown for simulations with varied cosmological 
parameters with that cosmologies' mass function using the multiplicity func- 
tion predicted from the central cosmology. Cosmological parameter are held 
constant, except triangles are for og = 0.8; squares are for Q m = 0.2; circles 
for h = 0.6 and a% = 0.8; diamonds are for Q„, = 0.2 and = 1 .2. 

A potential implication of this work is the calculation of 
the sensitivities of current and future cluster surveys to prop- 
erties of the primordial perturbations and the dark-energy- 
dependent growth of these perturbations. Some of the original 
work on quantifying parameter sensitivities of future cluster 
surveys considered t he bias associated with mass f unction un- 
certainties (iHolder et al II200U lHaiman et alJ|2001l) . as well 
as in more recent work (Battye & Weller 2003). However, the 
uncertainties in the predicted dark matt er halo ma ss function 
are often entirely ignored (Wang et al. 2004, 2005). 

To test the potential implications of dark matter halo mass 
function uncertainties on forecasts for parameter sen s itivity , 
we use the Fisher matrix approach of Holder et al. (2001) 
for sensitivities of the South Pole Telescope (SPT) Sunyaev- 
Zeldovich Effect cluster survey, with redshifts determined by 
the Dark Energy Survey (DES), and a minimum mass cut- 
off of M min = 3 x 10 14 h~ l M Q . Sky coverage is 4000 deg 2 
and cluster counts are binned in redshifts Az = 0.05 from 
z = to z = 2. We have also calculated and included the 
Fisher information from forecasts for the Planck mission 
parameter sensitivities from the cosmic microwave back- 
groun d (CMB) TT, TE and EE correlations (Eise nstein et al.l 
1998). In this test, we have not included scatter in the mass- 
temperature relati on (| Leyine e t alJl2002l) . nor uncertainty in 
this distribution ( Lima & Hu 2005|), bu t have also not in- 
clud ed mass-observable self-calibration ( Maiumdar & Mohr 
120031 which can be compensating effects. We find that the 
inferred uncertainties forecast for parameters best constrained 
by an SPT+DES cluster survey beyond the CMB constraints, 
namely the equation of state, fi m and neutrino mass, have 
forecast errors increased by, at most, approximately 1% in 



the error. Therefore, the mass function presented here has 
its uncertainties sufficiently well determined so that they are 
negligible compared to observational uncertainties. 

4. CONCLUSIONS 

We have measured the shape and quantified the uncertainty 
in the predicted theoretical mass function of dark matter ha- 
los in the ACDM cosmology with sixteen separate 1024 3 dark 
matter structure formation simulations. We found and cor- 
rected a systematic error in halo mass determina tion by the 
friend s-of- friends halo finder. T he canonical Shet h & Tormenl 
( 1999) and Jenkins et al. (2001) forms of the mass function 
are inconsistent with the ACDM mass function at ~10% level 
at intermediate masses, and >30% at the highest masses. 

Combining the work presented here with uncertainties in 
halo biasing will be necessary for impr oving work on quan- 
tified applications of the ha lo model llvan den Bosch et all 
120031 lAbazaiian et alJ 12005 ). The peak-background split 
model for the relation between halo bias and the mass function 
appears consistent, although this should be studied in further 
detail. 

Analyses of current and future cluster surveys' sensitivi- 
ties to the primordial perturbation spectrum and dark energy 
are not greatly dependent on inherent statistical uncertainties 
in the predicted dark matter halo mass function. The use of 
high-redshift cluster count surveys will require quantification 
of the predicted evolution of the mass function, as well as any 
cosmological dependence of the halo mass multiplicity func- 
tion. Simulation efforts may be tailored for a specific survey's 
requirements in this regard. In concert with observations, this 
will help elucidate the nature of dark energy, neutrino mass, 
cosmological structure and galaxy formation. 
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